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Abstract 

Using the quasiparticle self-consistent GW (QSGW) and local-density (LD) approximations, we calcu- 
late the g-dependent static dielectric function, and derive an effective 2D dielectric function corresponding 
to screening of point charges. In the g— >0 limit, the 2D function is found to scale approximately as the 
square root of the macroscopic dielectric function. Its value is ~4, a factor approximately 1.5 larger than 
predictions of Dirac model. Both kinds of dielectric functions depend strongly on q, in contrast with the 
Dirac model. The QSGW approximation is shown to describe QP levels very well, with small systematic 
errors analogous to bulk sp semiconductors. Local-field effects are rather more important in graphene 
than in bulk semiconductors. 
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Graphene is a first truly two-dimensional (2D) crystal, with unique electronic and structural 
properties (for review, see Refs. |l|-[5|). Screening of electron-electron and electron- impurity interac- 
tions in graphene is an important theoretical issue crucial for both many-body effects in electronic 
structure^ and for transport properties, especially, for electron scattering by charge impurities^i^. 
There are numerous works ''"^^ treating this issue within the two-band Dirac model. But the Dirac 
model does not take into account the many other bands involved, which can include van Hove 
singularities in electron density of states^ that may possibly be very essential, specifically for 
screeningi^ii^. Here we develop a definition for an effective 2D dielectric function in an ab initio 
context, and calculate it within the Quasiparticle self-consistent GW (QSGW) and local-density 
(LD) approximations. The former takes into account many-body effects beyond the density func- 
tional GGA or LDA schemes essential for correct description of excited states and thus screening 
effectsiS^. 

There are several GW calculations for graphene^Si"— , where G and the screened Coulomb in- 
teraction W are computed from the LDA. They all predict a notable (20-40%) increase of the 
Fermi velocity vp at the Dirac point K relative to the LDA(GGA) value, with vp between 1.1 
and 1.2-10^ m/s, in very good agreement with experiment^"-^. The dielectric function and optical 
condu ctivity as a function of frequency u for zero wave vector q=0 was also calculated in Refs. 
Here we focus on the static dielectric function {uj=0) as a function of q. This quantity 
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is relevant for calculations of resistivity via charge impurities^i^, as well as for the problem of 
supercritical Coulomb centers^'^^i'^^ and possible exciton instabilities^. 

The inverse dielectric function e~^(r, r', u) relates the change in total potential 6V to an external 
perturbing potential SV^^^ asr^^ 

6V{r,u) = j dr'e-\r,r',u)6V''''\r',u). (1) 

is obtained from a convolution of the polarization operator H and the bare Coulomb interaction 

V as 

e"^ = {l-vU)-\ 

In a system with translation symmetry, e~^, H, and v can be expanded in Bloch functions 
{Bf{r)h e.g. 

6-i(r,r',u;) = V ^Bf{r)ej}{ci,oo)Br{r') (2) 
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The most common choice of {Bf{r)} are plane waves, 

Bf{r) ^ B^{r) = exp(^(q + G) ■ r), (3) 

G being reciprocal lattice vectors. 

Quantities of interest are coarse-grained averages of e^^,{q^,uj). The "macroscopic" response 
to a plane wave perturbation is^ 

eM (q, w) = [eG=o,G'=o (q> ^)] ~^ (4) 

The matrix structure of e^^ with G 7^ G' reflects local field effects in terms of classical electro- 
dynamics. The quantity eAf(q) is commonly approximated by just e(q); that is, the Umklapp 
processes, or local field effects are neglected. This is not such a bad approximation in sp semicon- 
ductors but as we show here, it is a rather poor approximation in graphene. CmIq) corresponds to 
screening potential 6V^^^ with a single Fourier component q. Selecting G=G'=0 averages over 
the unit cell, restricting the spatial variation to the envelope exp(iq • r). While eA/ is a quantity 
of relevance to some experiments, perhaps the most relevant is screening of a point charge in the 
graphene sheet, which governs e.g., scattering from impurities. 

As graphene is a 2D system, we need to consider how the impurity potential v{q) = An/q'^ is 
screened in the sheet. The (statically) screened potential from a point charge at the origin may 
be written in cylindrical coordinates r={p,z,9) and (i={q,qz,(^q) as 

W{p,z) = — I dqqJo{qp)W'''{q,z) (5) 





W'^iq,z) = 4 r%e^^-^4%^ (6) 
Jo Qz + Q 

Thanks to graphene's hexagonal symmetry, does not depend on 6g for small q. W'^^{q,z) is 
the 2D (Hankel) transform of W{r), the analog of the 3D transform W{q)=e~^{q)v{q^). In the 
absence of screening = 1 and W'^^{q, z) reduces to the bare coulomb interaction z): 

v'^^iq, z)=A dqze''^^^—^ = —e'^^ (7) 
Jo Qz+Q Q 

An appropriate definition of an effective 2D dielectric function is then 

e2°(g,z) = ^2°(g,z)/iy2°(g,z) (8) 

Graphene wave functions have some extent in z which must be integrated over to obtain a scatter- 
ing matrix element. But the largest contribution originates from z=0, so W'^^{q, 0) is a reasonable 
estimate for the scattering potential. This is particularly so for small q of primary interest here. 



In practice we carry calculations in a periodic array of graphene sheets in the xy plane, spaced 
by a distance large enough that the sheets interact negligibly. To calculate eQLG'=o(^' 
adopt the all-electron, augmented wave implementation that was developed for the quasiparticle 
self-consistent GW (QSCVF) approximation, described in Ref.-^. It makes no pseudo- or shape- 
approximation to the potential, and does not use PWs (Eq. [3D for the product basis {-B}, but 
a mixed basis consisting of products of augmented functions in augmentation spheres, and plane 
waves in the interstitial region. The all-electron implementation enables us to properly treat 
core states. We calculate e^^(q, w) in the random phase approximation, using Bloch functions 
for eigenstatesi^. These are obtained from single-particle eigenfunctions \l/kn and eigenvalues 
in both the LDA and QSGVF approximations. In both cases the generalized LMTO method is 
usedSI'SS. 

QSGVF has been shown to be an excellent predictor of materials properties for many classes 
of compounds composed of elements throughout the Periodic Table, with unprecedented ability 
to consistently and reliably predict materials properties over a wide range of materials^'^^^— . 
Nevertheless there are small, systematic errors: in particular bandgaps in insulators such as GaAs, 
SrTiOa and NiO, are systematically overestimated. Its origin can be traced to a large extent to the 
RPA approximation to the polarizability, I]^^^=iGG. The RPA bubble diagrams omit electron- 
hole interactions in their intermediate states. Short-range attractive (electron-hole) interactions 
induce redshifts in Im e(q, u) at energies well above the fundamental bandgap; see e.g. Fig. 6 
in Ref.— ^. That ladder diagrams are sufficient to remedy most of the important errors in n^^"^ 
was demonstrated rather convincingly in CU2O, by Bruneval et al.— . Moreover Shishkin et al^ 
incorporated these ladder diagrams in an approximate way for several sp semiconductors, and 
established that they do in fact largely ameliorate the gap errors. Yang et al. investigated the 
effect of ladder diagrams in graphene and graphite, and showed that in a manner very analogous 
to ordinary semiconductors, these diagrams induce a redshift in the peak of Ime^^'^(q, near 
5 eV,-^ of ~0.6 eV. They found a strong similarity with conventional semiconductors, namely that 
the redshift from ladder diagrams approximately cancels the error in the LDA joint density of 
states. 

A redshift in the peak of Ime(ci;) increases the static dielectric constant 60^, as can be readily 
seen through the Kramers-Kronig relations. Remarkably, eoo as calculated by the RPA in QSCVF, 
is underestimated by a nearly universal factor of 0.8, for many kinds of insulators and semiconduc- 
tors'^°, including transition metal oxides such as NiO^-, Ce02, and sp semiconductors^. (This error 
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TABLE I: Energy gap Eq and valence bandwidth Tiy in diamond (eV); Fermi velocities vp in graphite 
and graphene (lO^m/sec). There is a significant renormalization of the bandgap from the electron-phonon 
interaction in diamond, estimated to be 370 meV— . Thus QSGW overestimates Eq by a slightly smaller 
amount than in other semiconductors, and the scaling of S as described in the text results in a slightly 
underestimated gap. The electron-phonon interaction also reduces the Fermi velocity in graphene, esti- 
mated to be 4 to 8% in an LDA-linear response calculation^^. The calculated Fermi velocities should be 
reduced by this much when comparing to experiment, vp calculated by QSGW is slightly overestimated, 
for much the same reason semiconductor gaps are overestimated, vp calculated from the scaled-S po- 
tential, is slightly larger than vp calculated LDA-based GW, i.e. G^^^W^^^^"^ , just as semiconductor 
bandgaps are slightly larger. When renormalized by the electron-phonon interaction, vp agrees very well 
with the measured valued. 

LDA QSGW scaled S Expt 

Eg, diamond 4.09 5.93 5.56 5.50 

Fi^, diamond 21.3 23.1 22.7 23.0 ± 0.2" 

Fi^, graphene 19.4 22.9 22.2 

?;f(H), graphite 0.77 0.99 0.94 0.91 ±0.15 

vp{K), graphene 0.82 1.29 1.20 1.1 



-^Ref. 
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is often approximately canceled in the LDA, fortuitously. As Yang et al. noted, the cancellation 
seems to apply to graphene in a manner similar to ordinary semiconductors.) Because e is system- 
atically underestimated, W = e~^v and S = —iGW are systematically overestimated; therefore 
QP excitation energies are also systematically overestimated. We have found that simply scaling 
by 0.8 (the nearly universal ratio e^*^"^ /e^^^) largely eliminates discrepancies between QSGVF 
and measured QP levels in a wide range of spd systems, including all zincblende semiconductors, 
and many other kinds of insulators. For graphene, we find that the QSGVF macroscopic (q'— J-O) 
dielectric constant was found to be 80% of the LDA one, consistent with the universal pattern 
in bulk insulators noted above. The many points of consistency with 3D behavior, both in the 
Q'&GW QP levels and the dielectric response suggest that QSGVF will exhibit the same reliable 
description of the 2D graphene system, with similar systematic errors. To confirm this, some band 
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FIG. 1: QS GW bands of graphene (dotted red lines), compared to LDA results (dashed blue lines) and 
QSGW results with S scaled by 0.8 (solid green lines) described in the text. The linear dispersion near 
K (or H, in graphite) is significantly larger in the QSGW case. Differences are quantified in Table HI The 
lowest lying unbound state can be seen as a parabolic band starting at T near 3.5 eV. It corresponds to 
the work function. LDA and QSGW work functions are very similar, consistent with the observation 
that LDA predicts work functions rather well in many systems. 



parameters for three pure (undoped) carbon compounds calculated by QSGW and QSGW with 
S scaled by 0.8 are shown in Table HI Scaling QSGW has a minor effect on the quasiparticle levels: 
e.g. it reduces vp by 7%. As Table [T] shows, vp falls in very close agreement with experiment when 
S is scaled and the electron-phonon interaction is taken into account, consistent with agreement in 
gaps in the bulk insulators. Even though the QSGW and LDA work functions are similar (Fig. [1]), 
the valence band is significantly widened relative to LDA,^ more so in graphene than in diamond. 

Careful checks for convergence were made in various parameters. To check supercell artifacts, 
a "small" 3D unit cell with the graphene planes repeated at a spacing equivalent to 4 atomic 
layers of graphite (25 a.u.) was compared against a "large" cell, with graphene planes spaced at 
8 layers. The bands from — oo to Ep+5 eV were found to be a very similar, with a slight increase 
in vp (1.23— )-1.29 -10^ m/s). k convergence in the construction of S was monitored by comparing 
QP levels generated on a 6x6x2 k mesh to a 9x9x2 mesh. QP levels were nearly identical: vp 
differed by <1% in the both the small and large 3D cells. 

eoo^(q, cj) must be integrated with a fine k mesh. To deal with the delicate q— )-0 limit, we 
calculated integrating on a standard k mesh including F, and an offset mesh (Eqns. 47 and 
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26l ). and averaged them. We present data for averaged 18x18x4 meshes. Calculations 
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without local fields were also performed for a pair of 24x24x4 meshes. e(q||, '?2=0, c<;=0) calculated 
by 18- and 24- (averaged) mesh integrations were essentially indistinguishable for g>0.1x27r/a, 
and differed by a few percent for g>0.02x27r/a. 
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FIG. 2: (Top) Static dielectric function eoo(^, ^2=0) along the (100) line in graphene, with local fields 
included and without, q is in units of 27r/a=2.56A^^. The q^O limit is delicate and there is some 
uncertainty in its value. Shown for comparison is the same function calculated from the LDA potential. 
In the ^^0 limit, eoo calculated by QSGW is -^O.S smaller than the LDA result, similar to the ratio found 
in bulk semiconductors. (Bottom) Effective layer dielectric function e^^(g, z=0) as defined by Eq. ([8]), 
with local fields, calculated within the QSGW and LD approximations. Local fields significantly reduce 
eoo- The LDA result for eoo(^=0.086,g'2=0) without local fields is ~4, which agrees with the w— t-O limit 
in Fig. 11 of Ref^- 



eoo{q, qz) was calculated on a grid of points {g, g^}; the ^2=0 case is shown in the first panel of 
Fig. O It was found that eoo is well parametrized (max error <0.1) by 

a?{q) + ql 



a\q) 



eoo(g,0)a2(g) + ql 

^2 



a^aiq 



aoq^ 



(9) 
(10) 



ai + q^ 

where ao=1.3 and 1.2 for QSGVF and LDA, respectively, and ai=1.6(27r/a)^. The approximate 
form for a in Eq. [10] is sufficient for any q where eoo differs significantly from unity. With Eq. ([9]) 
W'^^ can be integrated analytically. Taking the approximate expression for a^(g) we obtain 

7(7^ - 1) 
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7(ao - 1) + (7^ - ao)e{i-T)9^ 

where 7 = A/aoeoo(g, 0). 

Fig. |2] shows both kinds of dielectric functions, eu corresponding to the macroscopic polariz- 
ability, and the effective 2D static dielectric function e^°(g, 2;=0) calculated from Eq. Local 
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fields reduce the strength of the screening. The difference between LDA and QSGW results are 
modest; and as noted earlier, the LDA results are likely to be slightly better because they benefit 
from a fortuitous cancellation of errors. As g— j-O, 7 is significantly larger than oq and unity. Keep- 
ing only the leading order in 7, we obtain the surprising result that e^^(0,2;=0)~A/aoeoo(^5 ?2=0). 
e^^(0,z=0) is roughly a factor 1.5 times larger than the Dirac Hamiltonian result at zero doping. 
Such a model predicts e(g) ~ 2.4 independent of q, as shown by Andc^. We find e^°(g, 2=0)~3.5 
for but is a very strong function of q. 

Although virtual transitions involving Van Hove peaks of the density of states might strongly 
enhancei^ii^ were they sufficiently close to the Fermi level, apparently lie too far away in 
graphene. The case of small q (g~/ci?<10~^ A~^) is relevant for transport properties. In this region 
our first-principles calculations do not dramatically contradict predictions of the Dirac model. At 
the same time, for the problem of supercritical Coulomb centers and relativistic collapse (fall on 
the center distances of order of several inverse lattice constants are essential (this is the 
radius of screening cloud, according to renormalization group analysis^), which corresponds to 
larger q. For this region our results show that the Dirac model overestimates the screening. 
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